Thermomechanical properties of a single hexagonal boron nitride sheet 
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Using atomistic simulations we investigate the thermodynamical properties of a single atomic layer 
of hexagonal boron nitride (h-BN). The thermal induced ripples, heat capacity, and thermal lattice 
expansion of large scale h-BN sheets are determined and compared to those found for graphene 
(GE) for temperatures up to 1000 K. By analyzing the mean square height fluctuations {h'^) and 
the height-height correlation function H{q) we found that the h-BN sheet is a less stiff material as 
compared to graphene. The bending rigidity of h-BN: i) is about 16% smaller than the one of GE at 
room temperature (300 K), and ii) increases with temperature as in GE. The difference in stiffness 
between h-BN and GE results in unequal responses to external uniaxial and shear stress and different 
buckling transitions. In contrast to a GE sheet, the buckling transition of a h-BN sheet depends 
strongly on the direction of the applied compression. The molar heat capacity, thermal expansion 
coefficient and the Gruneisen parameter are estimated to be 25.2 JmoP^K"'^, 7.2xlO~®K~^ and 
0.89, respectively. 



I. INTRODUCTION 



A single layer of hexagonal boron-nitride (h-BN) is 
a wide gap insulator that is very promising for opto- 
electronic technologies [l|, 0], tunnel devices and field- 
effect transistors [SHa]. According to the well known 
Mermin- Wagner theorem Q, the stability of any two 
dimensional crystal is only possible in the presence of 
atomic corrugations which distort the perfect honey- 
comb lattice and allow the atoms to explore the out-of- 
plane direction. Experimental observations have found 
ripples in suspended sheets of GE [3, Q and atomistic 
simulations suggest that the strong bonds between the 
carbon atoms determine the features of these ripples [9j . 
Understanding the behavior of the ripples is important 
for many reasons [Ifll- They affect the electronic trans- 
port properties, e. g., in GE the ripples are believed to 
be one of the dominant scattering sources which limits 
the electron mobility [ll|, [l2] ■ 

h-BN ribbons doped by carbon has recently been pro- 
posed [3, [l^. In addition, BN based nanostructures 
are potential materials for thermal management applica- 
tions |14l4l3 | because of their high thermal conductivity 
and sensitivity to isotopic substitution and etc. There- 
fore, the knowledge of the shape and the temperature 
dependence of the intrinsic ripples is fundamental to de- 
vise novel nano-devices based on this material. 

Both GE and h-BN sheets have a honeycomb lat- 
tice structure, however the different electronic proper- 
ties and phonon band structure |20l - |22| results in un- 
equal morphologies and corrugations. Transmission elec- 
tron microscopy is widely used to resolve the individ- 




* corresponding author: ncckanial@srttu.edu 



Armchair 
O Boron O Nitrogen 

FIG. 1: (Color online) Schematic view of the single h-BN 
sheet. Smaller-yellow (bigger-blue) circles refer to the B (N) 
atoms. The rectangles indicate the atoms that are fixed dur- 
ing compression. Dashed (straight) lines correspond to arm- 
chair (zig-zag) uniaxial compression in the direction given by 
the arrows. Open arrows indicate the shear stress applied in 
the armchair direction. 



ual atoms in suspended h-BN sheets [23| where ripples 
inherently exist. First-principle calculations have been 
performed using small unit cells, periodically replicated, 
which are unable to model long wavelength corrugations 
which require thousands of atoms [2J| while the mechan- 
ical properties of a h-BN sheet can be estimated by us- 
ing DFT 123]. Conversely, atomistic simulations using 
molecular dynamics simulations (MD) enable to study 
thermo-mechanical properties directly on large scale sam- 
ples. The modified TersofF potential [26,] (TP) (parame- 
terized originally for carbon and silicon) with various set 
of parameters have shown to predict reasonable values for 



the thermo- mechanical properties of the h-BN sheet. In 
the pioneer work by Albe et al re-parameterized TP was 
used to study the impact of energetic boron and nitrogen 
atoms on a cubic-BN target 27]. Some other groups have 
also parameterized TP using various experimental data, 
e.g. Sekkal et al [2a] treated h-BN as a one-component 
system, using the same potential parameters for both 
boron (B) and nitrogen (N) (neglecting the B-N interac- 
tion) to investigate the structural properties. Matsunaga 
et al Q proposed the TP of B in order to simulate cu- 
bic boron carbonitrides which are cornpatible with the 
parameters of nitrogen fitted by Kroll [30| , and recently, 
Liao et al [31| and Sevik et al [32| reported TP param- 
eters that were used to study the mechanical properties 
and the thermal conductivity of a h-BN sheet, respec- 
tively. 

In this study, we investigate the thermal rippling be- 
havior of free standing monolayer h-BN by using state 
of the art molecular dynamics (MD) simulations. We 
adopted the TP potential re-parametrized by Sevik et 
al which has been shown to represent the experimen- 
tal structure and the phonon dispersion of the two- 
dimensional h-BN sheet. We found that h-BN is more 
corrugated and a less stiff material as compared to GE. 
The height-height correlations can be explained by the 
theory for continuum membranes [33]. In addition, we 
report results of both uniaxial and shear stress of a h-BN 
sheet and compare it with those found for GE. The buck- 
ling transition for compressed h-BN occurs earlier than 
for GE and the induced pattern of ripples when subjected 
to either uniaxial or shear stress shows significant differ- 
ences. 

This paper is organized as follows. In Sec. II, we in- 
troduce the atomistic model and the simulation method. 
Then, in Sec. Ill we analyze the behavior of the thermal 
ripples of a h-BN sheet. Here, we obtain the bending 
rigidity, the heat capacity, the thermal expansion coeffi- 
cients, and we study the buckling transition of the h-BN 
sheet when uniaxial external strain and shear stress are 
applied. All the results are compared with the ones ob- 
tained for a GE sheet. Finally, in Sec. IV, we present 
our conclusions. 



II. METHODS 

Classical atomistic molecular dynamics (MD) simula- 
tion is employed to simulate large scale samples of h-BN 
and GE at temperatures varying from 10 K up to 1000 K. 
We used a modified Tersoff potential (which is an ordi- 
nary defined potential in the LAMMPS package [33.l35|) 
according to the parameters proposed by Sevik et al for 
the h-BN sheet. All the parameters and a detailed de- 
scription of the potential energy is given in Ref . [32| . To 
simulate GE we have used the AIREBO potential [3g 
which is suitable for simulating hydrocarbons. We em- 
ployed NPT ensemble simulations with P=0 using the 
Nose-Hoover thermostat which enables us to allow the 



size of the system to equilibrate (i.e. the size of the 
system is not fixed). All the reported values have been 
computed averaging over 300 — 400 snapshots taken over 
uncorrelated samples. 

We start with a square shaped h-BN sheet with 
periodic boundary conditions and initial dimensions 
322A x32ll (315A xSlsA for GE) in the x and y direc- 
tion, respectively, which correspond to a total number 
of A^=37888 atoms, and which are sufficiently large in 
order to capture the long-wavelength regime. Periodic 
boundary conditions were adopted in both directions. 

To analyze the thermal ripples we computed the 
diffraction pattern which is typically studied in experi- 
ments to detect the size and shape of the corrugations [3| . 
We obtained also the mean square value of the out-of 
plane displacements JJh^) of the atoms and, by follow- 
ing previous works [33, |3g], the height- height correla- 
tion function {H(q)) which was determined by consid- 
ering an average of the height between the first neigh- 
bors of each atom. The latter was shown to follow a 
q^^ behavior that is expected from the theory of con- 
tinuum two-dimensional membranes at large values of q 
in the harmonic approximation (see below). To analyze 
the differences between the strain induced corrugations 
in the h-BN and the GE sheets we applied uniaxial and 
shear stress on both systems as is schematically shown 
in Fig. [T] In order to apply strain we set the force on 
the two ends equal to zero and move the end atoms with 
an infinitesimal compression step {dx = O.OlA) in the 
desired direction. After each compression step we wait 
for 2 ps to allow the system to relax and to ensure that 
the system is in thermal equilibrium [39|. Uniaxial com- 
pressive stress is applied along the zig-zag or armchair 
direction separately, and the shear stress is applied on 
the armchair edges with the opposite velocity for the top 
and bottom edges. The details of the used method of ap- 
plying the boundary stress can be found in our previous 
studies [H-Iill 

The TP function [H] used in the LAMMPS pack- 
age [SJ, I33 can be expressed as 
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where /c, /r and /a are cutoff functions, the repulsive 
pair term, and the attractive pair term, respectively. Ty 
and bij are respectively the distance from atom i to atom 
j and the bond order function. The use of TP disregard 
contributions coming from charge re-distribution which 
may occur in an ionic crystal. The inclusion of this effect 
in h-BN modifies the phonon spectrum significantly only 
for energies corresponding to the optical modes [42|, |43| . 
The large scale thermal ripples addressed here are gov- 
erned mainly by the transversal acoustic mode (TA), 
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FIG. 2: (Color online) (a) Contour plot of the heights for an arbitrary snapshot taken during the MD simulation at 300 K. 
(b) Corresponding simulated diffraction pattern for the sample shown in (a), (c) Zoom view of (b) around q=(0,0). 
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FIG. 3: (Color onhne) Variation of (a) < h^ > and (b) 
the bending rigidity versus temperature for both h-BN (open 
squares) and GE (solid circles), (c) Height-height correlation 
function H{q) of h-BN at two different temperatures as is 
indicated. The dashed lines are the harmonic q~* law. 



which accounts for out-of plane fluctuations, and it cou- 
ples with the in-plane modes which renormalizes the long 
wave-length behavior, e.g. the bendiiig rigidity k can be 
calculated directly from the TA mode l4J] . Therefore, the 
charge redistribution is expected not to affect the ther- 
mal fluctuations analyzed here and the use of the TP is 
justified ^^. 



III. RESULTS AND DISCUSSION 



A. Thermal excited ripples 



In Fig.[31Ja) we show a height contour plot of the atoms 
of the h-BN sheet for an arbitrary snapshot taken during 
the MD simulation at 300 K. The corresponding modeled 
diffraction pattern is shown in Fig.[2][b). This pattern is 
very similar to the one obtained for the GE sheet Q with 
the main difference in the distance between the Bragg 
points due to the unequal lattice constant of h-BN and 
GE. From the zoom plot around q=(0,0) (Fig.[2][c)) one 
observes the local structure of the central Bragg point 
for these intrinsic thermal ripples of h-BN. Notice that 
the lack of the presence of the q=(0,0) component is a 
consequence of the absence of a perfectly flat h-BN sheet. 



The signatures of the thermal induced ripples can also 
be seen in the mean square value of the out-of-plane 
fluctuations {h"^)- In Fig. ^a.) we show (/i^) as func- 
tion of temperature. In comparison with GE (included 
also in this figure for comparative purposes) we observe 
that (/i^) is always larger for the h-BN sheet. The 
weaker (stronger) the atomic B-N (C-C) bonds, the larger 
(smaller) the corrugations in h-BN (GE). Notice that the 
B-N bond is not pure covalent and it has an ionic charac- 
ter which is due to the different electronegativity between 
the two elements, i.e. the electrons are localized closer 
to the N atoms rather than the B atoms. Although this 
partially ionic structure of the h-BN layer increases the 
interlayer interaction resulting in a larger hardness of 3D 
bulk h-BN relative to graphite, it makes the single layer 
of BN less stiff than graphene. Moreover, the stacking 
of h-BN sheets in bulk h-BN is AAA stacking which is 
different for the ABC stacking in graphite j23i] . 



B. Bending rigidity K 



Accordingly to the two-dimensional theory of contin- 
uum membranes the height-height correlation function, 
in the harmonic approximation, where the out-of-plane 
and in-plane modes are decoupled, is expected to behave 
as 
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where k is the bending rigidity of the membrane, iV is 
the number of atoms of the sample, 5*0 the surface area 
per atom and fcs is the Boltzmann constant. In the large 
wavelength limit, i.e. for g — > 0, the height fluctuations 
are suppressed by anharmonic couplings between bend- 
ing and stretching modes giving rise to a renormalized q- 
dependent bending rigidity and hence Eq. ([3]) is no longer 
valid [Hill 0- In Fig. ^h) we show k for h-BN cal- 
culated from the harmonic part of H(q) between q=0.5 
A~^ and q=l A~^. In agreement with the larger values 
of (/i^), we observe that the h-BN membrane posses a 
smaller k as compared to GE and in the whole temper- 
ature range it is about 16 % smaller at room tempera- 
ture (300 K). Note that k for both h-BN and GE increase 
with temperature. In Fig. |3Kc) we show H{q) at 200 K 
and 1000 K were the harmonic behavior can be clearly 
observed (fitted with a dashed line) and, as expected, 
with increasing temperature H{q) is shifted to larger q. 
This figure also reveals that the ripples are not charac- 
terized by an unique wave-length and instead follows the 
behavior expected from continuum membrane theory. 

Before ending this section it is worthwhile to 
investigate an alternative method to estimate the 
bending rigidity (flexural rigidity). A common 
method for calculating the bending rigidity of 
BN-sheet is by performing simulations of BN- 
nanotubes as a function of its radius (R) of the 
curved tubes and then extrapolating the results 
to R^> oo. Hence, one can calculate the elastic en- 
ergy of the nanotube as a function of the inverse 
square of the radius, E ~ ^kR^^. This method was 
used in Ref. [25] and in our previous work [41] to 
calculate the zero temperature bending rigidity of 
GE and h-BN which were found to be 1.46 eV and 
1.29 eV, respectively. The result obtained with 
the Tersoff potential using Eq. (3) is less than the 
result of Ref. [25]. The bending rigidity of GE 
is larger than h-BN with about 0.15 eV in agree- 
ment with Ref. [25]. In order to have an indepen- 
dent check we estimated the bending rigidity of 
a BN sheet by performing several ground state 
simulations for (n,n) BN-nanotubes with n=5,.., 
20 using the Tersoff potential. Figure 4 shows the 
variation of the strain energy per atom as func- 
tion of the inverse square radius of the tube. Fit- 
ting a line to the data and dividing by the area of 
half of unit cell atom results in k = 0.86 eV. The 
result of Tersoff potential agrees well with our 
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FIG. 4: (Color online) Variation of strain energy versus 
inverse square of BN-nanotube radius using Tersoff 
potential. 
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FIG. 5: (Color online) Variation of (a) the total energy per 
atom and (b) the averaged bond length versus temperature 
for h-BN and GE (The error bar is less than 10~* eV/atom 
for the total energy). 



estimation for k using Eq. (3) (Ref. [25]). Thus 
we can conclude that result based on the Tersoff 
potential underestimate the bending rigidity as 
compared to the result from DFT. 



C. Heat capacity, thermal expansion and 
Gruneisen parameter 

The variation of the total energy per atom with tem- 
perature (> lOK) is shown in Fig. [5ja). The total en- 
ergy varies linearly with temperature and gives the corre- 
sponding lattice contribution to the molar heat capacity 
at constant pressure (the average size of the system after 



relaxation is taken constant) which for the h-BN sheet 
resuhs into 



Cp = 



dEr 
~dT 



= 25.2(3) J mor^iC-^ 



(4) 



which is comparable to the one for a GE sheet, i.e. 
24.5(9) J mol~^K~^ and close to the well known classical 
molar heat capacity, i.e., C = 3K ~ 24.94 J mol~^K^^, 
i.e., the Dulong-Petit value, where SR is the universal gas 
constant. Notice that the heat capacity is a tem- 
perature dependent parameter that will decrease 
for temperature below the Debye temperature. 
However, the classical MD simulation gives the 
correct high temperature asymptotic limit i.e., 
the Dulong-Petit value, but fails in the low and 
intermediate temperature range. 

In Fig. [5jb) we show the variation of the averaged B- 
N bond length with temperature. The linear behavior 
enables us to calculate the linear thermal expansion of 
the lattice parameter of a h-BN sheet as 



1 da 
a dT 



lBN = -^ = 7.2(1) xlO-^K-' 



(5) 



where a= 1.442 A is the zero temperature lattice parame- 
ter. The obtained 7 is 33% larger than the one measured 
for cubic BN, i.e. 4.8xlO~^K~^ [43| and is comparable 
to the one found for GE, i.e. jge = 10.0(7) x IQ-^ K-\ 
The latter (i.e. "(nE) is in good agreement with our pre- 
vious studies |39l.l40|. 

From the obtained values of 7 and C we can compute 
the Gruneisen parameter 



asN = -pr ^ 0.89, 
Cp 



(6) 



where B is the two dimensional bulk modulus and p 
is the mass density. Using Bh^BN=SeVk^'^ [Bqe = 
12.7 eVA-2 HI)), ph-BN = 2A.82/Sh-BN - 3.81 x 
lO^-^gm-^ [pQj^ = 7.6 X 10-4gm-2) for h-BN (GE), we 
obtain ubn = 0.89 (aGB=1.2). Note that the estimated 
value of aoE is better than the one found in Ref. |49|, i.e. 
-0.2, and is closer to the experimental result, i.e. 2.0 |50| . 



D. Buckling transition 

The different stiffness between the h-BN and GE sheets 
results in different buckling transitions. To study this, we 
compressed the systems uniaxially where we considered 
both zig-zag and armchair directions. This was done by 
fixing one row of atoms in each edge during the compres- 
sion steps as is indicated schematically in Fig. [TJ The 
compression rate was taken p = 0.5m/ s which is simi- 
lar to the one used in our previous study [40|, and small 
enough to guarantee that the system is in equilibrium 
during the whole compression process |5l| . The simula- 
tions were carried out at room temperature. Figure [BJ^a) 
shows the variation oi < h^ > versus strain, i.e. e = pt/l 
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FIG. 6: (Color online) Variation of < ft^ > versus external 
uniaxial strain in the (a) armchair and (b) zig-zag direction. 
Contour plot of the heights of compressed h-BN sheet samples 
subjected to a compressive strain for fixed < h^ >= 20 A^ 
for compression in (c) zig-zag direction of the h-BN sheet. 
Arrows indicate the stress direction. 



where t is the time (after starting the compression) and 
I is the initial box length in the direction of the compres- 
sion. The buckling transition for the h-BN sheet kept 
at room temperature is found to be 0.6 (0.1) which is 
smaller than the one for GE, i.e. 1.0 (1.2), for uniaxial 
compression along the zig-zag (armchair) direction 52]. 
Hence, when the compression is applied in the zig-zag 
direction, the h-BN sheet resists much more against the 
external uniaxial stress as compared to the case the stress 
is applied along the armchair direction. The smaller 
critical strain at which the buckling transition in 
the armchair direction of the h-BN sheet takes 
place indicates its more sensitive nature to exter- 
nal uniaxial stress [53]. Although the same argument 
holds for GE, the difference between the two directions 
is much smaller. Notice that DFT calculations re- 
sult in a deviation in bending rigidity (flexural 
rigidity) between ZZ BN-nanotube and AC BN- 
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FIG. 7: (Color online) Variation of buckling strain with the 
longitudinal length of the BN sheet which is compressed along 
the zig-zag direction. 



nanotube (ZZ becomes larger than AC) for radius 
smaller than ~ SA [25] while for larger radius they 
are the same. 

A contour plot of the buckled h-BN sheet with (ft,^)=20 
A^ and compressed in the zig-zag direction is shown in 
Fig. IHc). The obtained buckling transition for GE, i.e. 
0.8%, agrees very well with our previous studies [39| and 
is in the range of the experimental measured buckling 
transition for suspended GE, i.e. 0.5%-0.7% ^E^. 

It is also interesting to compare the buckling strain 
with those predicted by the theory of loaded beam. Euler 
theory for a two end-loaded beam having length L pre- 
dicts that efc oc L~^ [40, lla^. The first demonstration 
by MD of Euler buckling in nanostructures goes 
back to the pioneer work by Yakobson et al |.57J . 
We performed several simulations in order to find the 
length dependence of the buckling strain for BN sheets 
which are e.g. compressed along the zig-zag direction. 
Figure [7] shows the variation of ei, with the longitudinal 
length (L) of the zig-zag BN sheet. The solid line is a fit 
to L~^ and the symbols are the results of our MD sim- 
ulations which are qualitatively in agreement with the 
theory of loaded beam. 

Finally, we report the results for h-BN and GE sheets 



under shear stress. Here we applied the stress only in the 
armchair direction as is described schematically in Fig.[T] 
We found that, due to the larger stiffness, to reach the 
same value of {h"^) a larger shear stress has to be applied 
in GE as compared to h-BN. In Figs. |8l^a,b) one can ob- 
serve significant differences between the corrugated con- 
figurations of h-BN and GE sheets. These samples were 
subjected to a shear stress of e = 1.5%. While h-BN is 
highly sensitive to shear stress and deforms easily, GE 
can resist larger stress values due to its larger rigidity. 
IV. CONCLUSIONS 

The thermal properties of a boron nitride sheet were 
studied using large scale atomistic simulations. We 
showed that the scaling properties of a h-BN sheet fol- 
lows closely the results of membrane theory and hence 
the thermal excited ripples are not characterized by any 
particular wave-length. Using the harmonic part of the 
height-height correlation function we found an increas- 
ing bending rigidity with temperature which is smaller 
than the one of graphene. We found that the buckling 
transition for h-BN depends on the applied compression 
direction and is much smaller than the one of graphene. 
The obtained molar heat capacity agrees very well with 
the well-known Dulong-Petit number, 25.2 J mol~^K^^ 
and the thermal expansion coefficient was found to be 
positive and equal to 7.2 xlQ~^K~^. The Gruneisen 
parameter 0.89 is found to be smaller than the one for 
graphene, i.e. 1.2. We showed that the different stiffness 
between the GE and h-BN sheets leads to different pat- 
terns of deformations in the presence of either uniaxial 
or shear stress. 
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